local license =
[[----------------------------------------------------------------------
  **** BEGIN LICENSE BLOCK ****
	A simple flocking simulator written in LUA.

    Copyright (C) 2010
	Eric Fredericksen
	www.pttpsystems.com
	eric@pttpsystems.com

	This file is part of Flocker.

    Flocker is free software: you can redistribute it and/or modify
    it under the terms of the GNU Affero General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Affero General Public License for more details.

    You should have received a copy of the GNU Affero General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

  **** END LICENSE BLOCK ****
------------------------------------------------------------------------]]


-- Here is a simple function for dumping objects - a nice debugging tool
function edump(label, t, optionalFilter )
	if( not t ) then t = label; label = "dumping"; end
	if "table" == type(t) then
		for n, v in pairs( t ) do
			if not optionalFilter or optionalFilter(n, v) then
				print( tostring(label) .. ": ", n, v)
			end
		end
	else
		print(tostring(label)..": ", type(t), tostring(t) )
	end
end

function CopyArray( thisArray )
	local newArray = {}
	for k,v in ipairs(thisArray) do newArray[k] = v end
	return newArray
end

function CopyTable(thisTable)
  local newTable = {}
  for k,v in pairs(thisTable) do newTable[k] = v end
  return newTable
end

function AddToTable(thisTable, ...)
  for k,v in pairs(arg) do table.insert(thisTable, v) end
  return thisTable
end

--[[ ================= table serializer and deserializer ========================
	THIS IS NOT GENERAL PURPOSE. I use it for storing blobs in the database.
	It works for tables without circular references or shared references.
		Circular references will blow it up.
		Shared references will create duplicates upon reconstition.
		]]
function decodeLSONTable( _string )
	local reconstitutorFn = assert( loadstring( _string ) )
	return reconstitutorFn()
end

function encodeLSONTable( _table, recurse )


	local lson = "if true then local theTable = { " -- non recurse condition

	if recurse then lson = "{" end

	local thefirst = true
	for k, item in pairs( _table ) do

		-- only process the non functions
		if "function" ~= type(item) then

			if not thefirst then lson = lson .. ", "; else thefirst = false end

			if "number" == type(k) then

				if k == math.floor(k) then
					lson = lson .. string.format("[%d] = ", k)
				else
					lson = lson .. string.format("[%f] = ", k)
				end
			else
				lson = lson .. string.format("[%q] = ", tostring(k) )
			end

			if "table" == type( item ) then
				-- recure the table
				lson = lson .. encodeLSONTable(item, true)

			elseif "string" == type( item ) then
					lson = lson .. string.format("%q", item)
			else
				lson = lson .. tostring(item)
			end
		end

	end

	if recurse then
		lson = lson .. " }"
	else
		lson = lson .. " }; return theTable; end"
	end

	return lson

end

-- this needs to work for openGL, where pitch is rotation about X and yaw is rotation about Y
function PitchAndYawToCartesian( pitch, yaw, length, degrees)

	length = length or 1.0

	-- convert to radians
	if degrees then
		pitch 	= pitch * (math.pi/180)
		yaw 	= yaw   * (math.pi/180)
	end

	-- start with vector on the X-axis
	local x,y,z = length,0,0

	-- correct pitch sign from atan2()
	pitch = -pitch

	-- pitch about Y
	local cosY, sinY = math.cos(pitch),   math.sin(pitch)
	x, y, z = cosY * x + sinY * z, 		y, 							-sinY * x + cosY * z

	-- yaw about Z
	local cosR, sinR = math.cos(yaw), math.sin(yaw)
	x, y, z = cosR * x - sinR * y, 		sinR * x + cosR * y, 		z

	return x,y,z

end


function OGLCartesianToPitchAndYawOGL( ex, wye, zee, degrees )

	local length = (ex^2 + wye^2 + zee^2)^0.5
	-- normalize vector so we can compute angles
	ex, wye, zee = ex/length, wye/length, zee/length

	-- for computing angle with y-axis
	local el = (zee^2 + ex^2)^0.5
	local pitch = math.atan2( wye, el ) -- rotation around x-axis in y-z plane

	local yaw   = math.atan2( ex, zee ) -- rotation around y-axis in x-z plane

	if degrees then
		yaw 	= yaw   * (180/math.pi)
		pitch 	= pitch * (180/math.pi)
	end
	-- negate pitch cuz of our coord system
	return -pitch, yaw, length
end


function CartesianToPitchAndYaw( ex, wye, zee, degrees )

	local length = (ex^2 + wye^2 + zee^2)^0.5
	-- normalize vector so we can compute angles
	ex, wye, zee = ex/length, wye/length, zee/length
	-- for computing angle with z-axis
	local el = (ex^2 + wye^2)^0.5
	local pitch = math.atan2( zee, el  ) -- rotation around axis in x-y plane
	local yaw   = math.atan2( wye, ex  ) -- rotation around z

	if degrees then
		yaw 	= yaw   * (180/math.pi)
		pitch 	= pitch * (180/math.pi)
	end

	return pitch, yaw, length
end

-- unit testing for yaw-pitch <-> cartesian conversions
if false then
	local function EqualUpToPrecision( places, number1, number2 )
		local fmt = string.format("%%.%df",places)
		number1 = tonumber( string.format(fmt, number1) )
		number2 = tonumber( string.format(fmt, number2) )
		return number1 == number2
	end
	local tests =
	{
		 { 0,  0,  1 }
		,{ 0,  0, -1 }
		,{ 0,  1,  0 }
		,{ 0, -1,  0 }
		,{ 1,  0,  0 }
		,{-1,  0,  0 }
		,{ 1,  1,  math.sqrt(2) }
		,{-1, -1, -math.sqrt(2) }
		,{ 1,  0,  1 }
		,{-1,  0, -1 }
		,{ 1,  1,  0 }
		,{-1, -1,  0 }
	}

	for __,test in pairs(tests) do

	local p,y,l = CartesianToPitchAndYaw( test[1], test[2], test[3], true)
	local ex,wye,zee = PitchAndYawToCartesian( p,y,l, true)

	if false == (EqualUpToPrecision(15, test[1], ex)
					and EqualUpToPrecision(10, test[2], wye)
					and EqualUpToPrecision(10, test[3], zee)
					)
		then
			print( string.format("test:%10.4f %10.4f %10.4f -> %10.4f %10.4f %10.4f", test[1], test[2], test[3], p,y,l)  )
			print( string.format("pyl: %10.4f %10.4f %10.4f <- %10.4f %10.4f %10.4f", ex,wye,zee, p,y,l)   )
			print()
		end

	end

end


function CartesianToSpherical( ex, wye, zee, degrees )

	local x,y,z = ex, wye, zee

	if "table" == type(ex) then x,y,z = ex[1], ex[2], ex[3] end

	local r 	= (x^2 + y^2 + z^2)^0.5
	if r == 0 then return 0,0,0 end
	-- convert to degrees
	local theta	= math.acos(  z/r )
	local phi	= math.atan2( y, x )

	if degrees then
		theta	= (180/math.pi) * theta
		phi		= (180/math.pi) * phi
	end

	return r, theta, phi
end

function SphericalToCartesian( r, phi, theta, degrees )
	-- convert to radians
	if degrees then
		phi 	= phi   * (math.pi/180)
		theta 	= theta * (math.pi/180)
	end

	local sinPhi	= math.sin( phi )
	local x			= r * sinPhi * math.cos( theta )
	local y 		= r * sinPhi * math.sin( theta )
	local z			= r * math.cos( phi )
	return x,y,z
end

function SphericalToCartesianOddBall( r, theta, phi, degrees )
	-- convert to radians
	if degrees then
		phi 	= phi   * (math.pi/180)
		theta 	= theta * (math.pi/180)
	end

	local sinTheta	= math.sin( theta )
	local z, y, x	= 	r * math.cos( theta )
					,	r * sinTheta * math.sin( phi )
					,	r * sinTheta * math.cos( phi )
	return x,y,z
end

function RaisedCosine( now, period )
	-- smooth transition from 0.0 to 1.0 over period.
	-- smooth decelaration from 1.0 to 0.0 from period to 2*period
	local radians = -math.pi/2.0 + math.pi*now/period
	return  0.5 + 0.5*math.sin( radians )
end

-- http://en.wikipedia.org/wiki/Filter_%28higher-order_function%29
function FilterList( list, filterFn )

	local filtered = {}
	local count = 0
	for __, listItem in pairs( list ) do
		if not filterFn or filterFn( listItem ) then
			table.insert( filtered, listItem )
			count = count + 1
		end
	end
	return filtered, count
end



-- utilities for vector math
vec3 =
{
	  RADPERDEG = math.pi/180.0
	, DEGPERRAD = 180.0/math.pi

	, new   = function( x, y, z )
			return { x  or 0, y or 1, z or 0 }
		end
	, newV  = function( v ) return { v[1]  or 0, v[2] or 1, v[3] or 0 } end

	, copyTo = function( v1, v2 )
			v1[1], v1[2], v1[3] = v2[1], v2[2], v2[3]
			return v1
		end

	, equal  = function( v1, v2 )
			return v1[1] == v2[1] and v1[2] == v2[2] and v1[3] == v2[3]
		end

	, set = function( v, x, y, z )
			v[1], v[2], v[3] = x,y,z
			return v
		end

	, isNull    = function( v ) return 0==v[1] and 0==v[2]  and 0==v[3] end
	, isNotNull = function( v ) return 0~=v[1] or  0~=v[2]  or  0~=v[3] end

	, scaleIt = function( v, k ) return vec3.scale(v,k,true) end
	, scale = function( v, k, inplace )
			local x,y,z = v[1]*k, v[2]*k, v[3]*k
			if( inplace ) then
				v[1], v[2], v[3] = x,y,z
				return v
			end
			return {x,y,z}
		end


	, wrapIt = function( v, mx, my, mz ) return vec3.wrap( v, mx, my, mz, true ) end
	, wrap = function( v, mx, my, mz, inplace )
			mz = mz or 1
			local x,y,z = (v[1]+mx)%mx, (v[2]+my)%my, (v[3]+mz)%mz
			if( inplace ) then
				v[1], v[2], v[3] = x,y,z
				return v
			end
			return {x,y,z}
		end

	, addTo = function( v1, v2 ) return vec3.add(v1, v2, true ) end
	, add = function( v1, v2, inplace)
			local x,y,z = v1[1] + v2[1], v1[2] + v2[2], v1[3] + v2[3]
			if( inplace ) then
				v1[1], v1[2], v1[3] = x,y,z
				return v1
			end
			-- else
			return {x,y,z}
		end

	, subFrom = function( v1, v2 ) return vec3.sub(v1, v2, true ) end
	, sub = function( v1, v2, inplace )
			local x,y,z = v1[1] - v2[1], v1[2] - v2[2], v1[3] - v2[3]
			if( inplace ) then
				v1[1], v1[2], v1[3] = x,y,z
				return v1
			end
			-- else
			return {x,y,z}
		end

	, alphaBlendIt = function(v1, v2, alpha, inplace) return vec3.alphaBlend(v1,v2,alpha,true) end
	, alphaBlend   = function(v1, v2, alpha, inplace)
			local x = v1[1]*alpha + (1-alpha)*v2[1]
			local y = v1[2]*alpha + (1-alpha)*v2[2]
			local z = v1[3]*alpha + (1-alpha)*v2[3]
			if( inplace ) then
				v1[1], v1[2], v1[3] = x,y,z
				return v1
			end
			-- else
			return {x,y,z}
		end

	, dot = function( v1, v2 )
			return v1[1] * v2[1] + v1[2] * v2[2] + v1[3] * v2[3]
		end

	, distanceSquared = function( v1, v2 )
			return (v1[1] - v2[1])^2 + (v1[2] - v2[2])^2 + (v1[3] + v2[3])^2
		end

	, distance = function( v1, v2 )
			return ((v1[1] - v2[1])^2 + (v1[2] - v2[2])^2+ (v1[3] - v2[3])^2)^0.5
		end

	, length = function( v )
			return ( v[1]*v[1] + v[2]*v[2] + v[3]*v[3] )^0.5
		end

	, normalizeIt = function( v ) return vec3.normalize(v, true) end
	, normalize = function( v, inplace )
			local length = vec3.length( v )
			if( length < 1e-4 ) then

				-- arbitrary vector in case of zero length - pick one
				v[1] = 1/math.sqrt(2)
				v[2] = v[1]
				v[3] = 0
				length = 1
			end
			if( inplace ) then
				v[1], v[2], v[3] = v[1]/length, v[2]/length, v[3]/length
				return v
			end
			return { v[1]/length, v[2]/length, v[3]/length }
		end

	, rotDegIt = function( v, roll, pitch, yaw )
			return vec3.rot( v, roll, pitch, yaw, true, true )
		end

	, rotDeg = function( v, roll, pitch, yaw, inplace )
			return vec3.rot( v, roll, pitch, yaw, false, true )
		end

	, rotRadIt = function( v, roll, pitch, yaw )
			return vec3.rot( v, roll, pitch, yaw, true, false )
		end

	, rotRad = function( v, roll, pitch, yaw, inplace )
			return vec3.rot( v, roll, pitch, yaw, false, false )
		end

	-- rotation about the 3 axes
	, rot = function( v, roll, pitch, yaw, inplace, degrees )

			roll 	= roll 	and degrees and roll	* vec3.RADPERDEG
			pitch 	= pitch and degrees and pitch	* vec3.RADPERDEG
			yaw 	= yaw 	and degrees and yaw	* vec3.RADPERDEG

			local x,y,z = v[1], v[2], v[3]


			if pitch then	-- rotation about X
				local cosP, sinP = math.cos(pitch), math.sin(pitch)
				x, y, z = x, 						cosP * y - sinP * z,  		sinP * y + cosP * z
			end

			if yaw then	-- rotation about (new) Y
				local cosY, sinY = math.cos(yaw), math.sin(yaw)
				x, y, z = cosY * x + sinY * z, 		y, 							-sinY * x + cosY * z
			end

			if roll then	-- rotation about (new) Z
				local cosR, sinR = math.cos(roll), math.sin(roll)
				x, y, z = cosR * x - sinR * y, 		sinR * x + cosR * y, 		z
			end

			if( inplace ) then
				v[1] = x
				v[2] = y
				v[3] = z
				return v
			end
			return {x,y,z}
		end

			-- rotation about the z axis
	, rotOrig = function( v, phi, theta, inplace, isRads )

			if not isRads  then
				phi = phi * vec3.RADPERDEG
			end
			--[[ in the future, we might convert to spherical coords, adjust, then back again.]]
			-- if we are doing 3D then rescale the modification
			if theta then phi = phi/3 end
			local cosT, sinT = math.cos(phi), math.sin(phi)
			local x,y,z = v[1], v[2], v[3]
			-- rotation about Z
			x, y, z = cosT*x - sinT*y, 		sinT*x + cosT*y, 		z
			if theta then
				-- rotation about Y
				x, y, z = cosT*x + sinT*z, 		y, 						-sinT*x + cosT*z
				-- rotation about X
				x, y, z = x, 					cosT*y - sinT*z,  		sinT*y + cosT*z
			end

			if( inplace ) then
				v[1] = x
				v[2] = y
				v[3] = z
				return v
			end
			return {x,y,z}
		end



	, jitterPos		= function( v, from, to)
			if not (from and to) then from = -1 ; to =  1 end
			v[1] = v[1]+math.random(from,to)
			v[2] = v[2]+math.random(from,to)
			return v
		end

	, vecFromDeg	= function( phi, theta ) return vec3.vecFromTheta( phi, theta, false ) end
	, vecFromRad	= function( phi, theta ) return vec3.vecFromTheta( phi, theta, true  ) end
	, vecFromTheta	= function( phi, theta, isRads )

			if not isRads then phi   = phi   * vec3.RADPERDEG end

			if not theta then
				return { math.cos(phi), math.sin(phi), 0 }
			end

			if not isRads then theta = theta * vec3.RADPERDEG end

			return {
				  math.sin(theta)*math.cos(phi)
				, math.sin(theta)*math.sin(phi)
				, math.cos(theta)
				}
		end

}
